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We provide a general methodology to directly determine the validity of the meanfield Bogoliubov-de Gennes 
equation. In particular we apply this methodology to the case of two component interacting ultracold Fermi 
gases. As an example, we consider the case of population imbalance, between the two components, in 
the strongly attractive interacting regime, where meanfield results predict Fulde-Ferrell-Larkin-Ovchinnikov 
(FFLO) states. For these states we find at finite temperatures that the assumptions used to derive the Bogoliubov- 
de Gennes equation are invalid. 
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The onset of experimentally realized superfluidity of 
Fermionic gases by the MIT group in 2004 [1] has seen a 
trend towards theoretical studies of ultracold gases. With the 
use of a wide range of techniques at experimentalists finger- 
tips, it is possible to tune ultracold gases to sit in specific po- 
tentials, have precise particle numbers and create mismatched 
Fermi surfaces giving rise to many new concepts to study and 
analogies to solids [2]. An important advancement was the 
use of Feshbach resonances in controlling the s-wave scatter- 
ing length, flj; this was the defining step in realizing super- 
fluidity in two component fermionic systems. Applying an 
external magnetic field enables one to control the relative en- 
ergy between the two channels of the scattering process: the 
open (scattering) and closed (bound) channels. This gives ex- 
perimentalists the ability to form a smooth transition from the 
deep attractive to the strongly repulsed; the deep BCS side to 
tightly bound Bose-Einstein Condensation (BEC) interaction 
[3-6]. 

Almost as soon as superfluidity was achieved in Fermi 
gases the question of how spin population imbalance would 
effect the system was asked. When a gas has a polarization 
P - [N^ - Ni]/N greater than zero where is the total par- 
ticle number, a difference in the chemical potential for each 
spin state is enforced by the system, which can be thought of 
equivalently as applying an external magnetic field. Build- 
ing upon the work down by Bardeen, Cooper and Schrieffer 
[7, 8], the first analogous work done in superconductors was 
by Clogston and Chandresekaer [9, 10] where they proposed 
a critical cut-off' for the superconductivity in the presence of a 
magnetic field. 

An original idea for imbalance in superconductors was pro- 
posed in the 1960's by Fulde and Ferrel [11] and separately 
Larkin and Ovchinnikov [ 12], they found a spin density imbal- 
ance in a system would induce non-zero momentum Cooper 
pairs leading to oscillations in the pairing field, the Fulde, Fer- 
rel, Larkin and Ovchinnikov (FFLO) state. The experimen- 
tal realization of population imbalance in Fermi gases [13- 
16] in 2006 caused a flurry of theoretical investigations into 
the FFLO state. These have primarily focused on meanfield 
Bogoliubov-de Gennes approaches [17-22], which have been 
complemented by beyond meanfield Density Functional The- 
ory [23], variational techniques [24] and the T-matrix scheme 
[25] in the BEC-BCS crossover regime and ID exact and 
Quantum Monte Carlo techniques [26, 27]. This theoretical 



endeavor is complemented by experimental research [28-30] 
to observe the FFLO state. 

In this work we derive a general methodology for determin- 
ing the validity of the meanfield solutions of the Bogoliubov- 
de Gennes equation and implement this for the specific case of 
the well studied imbalanced ultracold Fermi gas in the strong 
interacting regime. As such, the paper is set out in the fol- 
lowing manner, firstly we describe the Bogoliubov-de Gennes 
meanfield theory and its regimes of validity. Then we present 
results for a balanced gas and an imbalanced gas in a spheri- 
cal trap, demonstrating, qualitatively, the regimes under which 
the meanfield solutions are valid. Finally, we comment on the 
validity of previous meanfield analysis of FFLO states in im- 
balanced ultracold Fermi gases. 

Our starting point for the study of balanced and imbalanced 
ultracold Fermi gases is the following microscopic Hamilto- 
nian 

cr ^ 

+ J dVdV'4']:(r')4'|(r)t/(r-r')4'i(r)4't(r'), (1) 

where 4'i-(r), 4'o-(r) are the real space annihilation and cre- 
ation operators for an atom with spin cr at position r. 'Ha = 
{-ft^V^ /Inicr - l-ia- + V'o-(r)) is the single particle Hamiltonian 
for a atom of mass nia- residing in a potential Va-{r), where 
jUo- is the chemical potential for component cr. U{r - r') - 
U6(r-r') is the contact interaction with U = AnTi^a Jnir, 
where niy is the reduced mass. 

Applying the meanfield approximation to Eq (1) and defin- 
ing the pairing field and spin density, A(r) = f/{*l'|(r)4'-f(r)) 
and no-(r) = (4'^(r)4'o-(r)) respectively, we can write the 
meanfield Hamiltonian as 

Hmf -^c + I d3r{24'^(r)'H^4'^(r) 

+ ^ U%{r)nAr)^Ar) + 4':f(r)A(r)4'|(r) + H.c], 

CT 

(2) 

where Qc is the ground state energy. 

It is well known that we can write Eq (2) in a diagonalized 
form, Hmf - Tjko- Ekyl^yf.^ by a Bogoliubov-Valatin trans- 
formation and solving the Bogoliubov-de Gennes equations 
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A(r) W «t(r) \ / M,(r) 
A*(r) -'Wi 11 v,(r) ) ""M v,(r) 



(3) 



where the wavefunctions for a particle with spin cr, position r 
and time t are given by 



k 



and are fermionic quasiparticle creation and annihi- 
lation operators and MA(r) and VA(r) are the Bogoliubov quasi- 
particle wave functions, with corresponding eigenenergies Ei^. 
We have introduced for completeness a Hartree term to the 
single particle Hamiltonian "Her = -fi^V- /2m,j- -/y,^ + Vtj-Cr)) + 
f/no-(r)«o-(r). The self-consistent pairing field A(r) and spin 
densities n,j{r) in terms of the quasiparticle wave functions 
are given by 



nT(r)=2|M,(r)|V(£A), 

k 

A(r)=f7(r)_2M,(r)v*(r)/(£,), 



(6) 
(7) 
(8) 



where the Fermi distribution is f(Ek) = [exp(-Ek/kBT) + 1]"' 
and U (r) has been introduced to renormalize the ultraviolet di- 
vergences, discussed in detail later These equations form the 
basis of most meanfield calculations of interacting ultracold 
Fermi gases, where the calculations are iterated to achieve 
self-consistency. 

To move from Eqs. ( 1) to (2) a meanfield approximation has 
been employed to the interacting part of the Hamiltonian and 
expanded about the small order parameters 

6(r, t) = U{r) [❖tCr, t)^^{r, t) - {^^{r, t)%{r, t))] , (9) 
SnAr, t) = %{r, t)^Ar, - {%{r, f)*.(r, O) , (10) 

where we have assumed the second order terms of 
(S"'"(r, f)S(r, f')) and {6nl-(r,t)6na-(r,t')} go to zero. We can 
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however calculate these terms retrospectively by defining the 
second order fluctuations in the pairing and density fields, for 
example in the pairing field 

<5'(r, t)6(r, t')) = [U{r)f [(4'|(r, f)*|(r, f)*i(r, t')%{r, t')) 

- {^{r, t)¥^{r, t)) (4'i(r, t')%{r, t')) 

(11) 

In an analogous method for current/charge fluctuations in 
mesoscopic conductors [32-34], the anomalous fluctuations 
are determined through the quantum statistical expectation 
value of the relation [35] 



<A(r, w)A(r, a>')) 



<A(r,w))<A(r,w')> 



(12) 



where A(r, w) = *l'|(r, w)*I'|(r, w) is the Fourier trans- 
form of the product of time space wavefunctions A(r, f) = 
4'|(r, t)^'^{r, t) and similarly to Eq. (9) we write 



6{r, CO) = U{r) A(r, w) - (A(r, w)) 



(13) 



Now we calculate the second order fluctuations in frequency 
space. Taking the Fourier transform of the product of wave- 
functions using the Bogoliubov quasiparticle operators we 
have 

A(r, w) =27r ^luikU^k+tiwJiky^fk+tm + v*k'^^k+nujj\{i'n+i>to 
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yiky\k-n., - ''*^k''lk-na>y]kylk-t,o\^ (i^) 



"ik^'{k-tm 



where the position r is implied; we then get a similar expres- 
sion for A*(r, w). Using the identity (Eq. (1.15))from [32], 

{yLynnylkysi) - {ylmypnXy'ykysi) 

= 6{E,„ - E„)6{Ek - E,)5,s5pyfa(E„)[l - fy{E,„)l (15) 

with Eq. (14) in Eq. (12) gives delta functions of the form 
6{Eif - Ej - TiLj') and 6{Ek + - Ej), allowing us to perform 
one summation over the energies and set w - -co' , giving 



{6\r)5{r))^ ^YJ[''*k-t^u?'^kV^k'^^k-nJ'^E^k-nJf{E'^k) + 2u*^_J^^^^ - f{E^k)) 

k 

«n+A.^/(%+^-)(l -/(%))]' 



}k+nw^'}k^*k^^^k+hJ^E^k+twj)f{E^k) + 2m 



^k+hiL}']k^'^k 



(16) 



and in the limit w — > we have the result. We wish to look at the zero frequency limit as the pairing field, 

Eq. (8) (and densities Eq. (6, 7)) are calculated in the zero 

{5\r)m) = Stt^ [U(r)f ^ \uu{rf\v,{rff{E,)f{-E,). 

k 

(17) 
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frequency limit. We can calculate in a similar fashion to above 
and find the anomalous fluctuations in the densities of each 
spin state 

<5«T(r)<5«T(r)) = ^7? Yj \iik{r)tf{Ek)f(-Ek), (18) 

k 

{6n^{r)5n^{r)) = Stt^ ^ \vu{r)t f{E,)f{-E,). (19) 

k 

In this paper we consider a spherically symmetric trap, i.e. 
y(r) = ma-apT^ 12, with the atomic mass of each of the com- 
ponents being equal, i.e. mn - m-\ - mi. For such a case 
it is natural to expand the wave functions 4'o-(r) in the basis 
of the spherical harmonic oscillator eigenstates. In the usual 
quantum numbers, following Refs [17, 20], [k] - {n,m,l] we 
have 

^a{r) = ^ gmnl{r)a„mla-, (20) 
nml 

where anmia- are the fermionic operators, gnmiir) = 

R.,ir,- V2(».„-)"'^;;;^^"?'C'"(?^)=-* (2.) 

are the normalized basis states, with t'^^^^ir^) being a La- 
guerre Polynomial and f = r/o/,,,, with = ^^jWJnv^iu be- 
ing the harmonic oscillator length. As w is a good quantum 
number we can calculate the angular integrations by introduc- 
ing a (2/ -H 1) degeneracy to the / states. The Z-states of the 
Hamiltonian can be diagonalized separately due to the axial 
symmetry, greatly reducing the computational complexity of 
the problem. In order to diagonalize the Hamiltonian and cal- 
culate any quantities we must truncate the quasiparticle en- 
ergies to Eni < Ec, where the choice of E^ is large enough 
to not change the final solution. Writing the order parameter 
and densities in this basis and introducing the I specific cut-off 
Ni ^{E^- I -3/2)/2 we obtain 

Nl 21+1 

A(r) = U{r) ^ -—R„,(r)R„,,(r){al,^al,^), (22) 

nn'l 

^21+1 .,. 
nAr) = 2j Rni(r)Rn'i{r){ali^a,fio-}. (23) 

nn'l 

These self consistent equations are solved by iteration for 
fixed particle numbers A^^- (majority) and Ni (minority) by 
adjusting the chemical potentials fia-- In each iteration the 
Hamiltonian is diagonalized and the resultant eigenvectors 
and eigenenergies are used to calculate the pairing and num- 
ber densities, Eqs. (22) and (23), until the relative change in 
the particle number and pairing field is < 10"^ between con- 
secutive iterations. 

We turn our attention now to the problem of ultraviolet di- 
vergences in the pairing function and Hartree terms. The di- 
vergence of the Hartree terms, Un^r) is clear as we go to the 



unitary regime a, — » oo, then J7 — > oo making the Hartree 
terms infinitely small. We can overcome this divergence by 
eliminating the Hartree terms from the free Hamiltonian in the 
unitary regime, and in our calculations we have not included 
them [36] . The divergence of the pairing field due to the nature 
of the contact interaction is easily solvable by regularizing the 
equation using the pseudo potential prescription [37], with the 
inverse of the regularized contact interaction given by 

1 I mu (k^(r) ^c(r) + ^F(r) \ 
^^^"^i^'" ^e(r)-Mr) -^-^^t ^''^ 

where the local Fermi momentum is given by ju = 
h^kF(rf/2m^l + V(r), E, = ft^k^(r)/2m-fi + V(r) and // = 
[jU| + iii]/2. Following from Ref. [18] we have used a hy- 
brid scheme for calculating the gap equation, i.e. A(r) = 
^BdG(j.) + ALDA(r), where A'^'*^(r) is given by Eq. (22) and 
^LDA^j.-j calculated from the semi-classical solution to the 
pairing equation obtained from the local density approxima- 
tion (LDA) for contributions above the cutoff E^. Specifically, 

^^], (25) 

w ith - E{r,k) + 6n, 6fi = - ni]/2, E(r,k) = 
V(e(r, k) - fif- + |A(r)p and e(r, k) = k^n^ /2mji + V(r). Using 
the hybrid scheme we found the solutions converge for a lower 
cut-off, speeding up our calculations. We can also include the 
relative fluctuations from the higher order corrections to the 
LDA approximation, however, any additional relative fluctua- 
tions will be positive definite, making ((5* (r)5(r)) larger Thus 
Eq. (17) provides a lower bound for (5' (r)(5(r)). 

Now that we have everything ready to begin solving the 
self consistent equations we can calculate our anomalous sta- 
tistical averages, Eqs (17) and (18), in the harmonic oscil- 
lator basis. Writing the quasiparticle eigenvectors in terms 
of the eigenstates a„i{r) = 2„ Rniir)W'j^^ and /3„i{r) = 
2„ Rni{r)W'j^^j^^ we have 

(dHmr)} = [U(r)fYj27r(2l+ l)a,ni(r)%i(rff(Eji)f(-Eji) 
jl 

(26) 

{6ni(r)6ni(r)} = ^ 2^(2/ + l)amiir)''f(Ej,)f(-Ej,) (27) 
jl 

{6n^ir)6n^(r)} = ^ 2n(2l + m„,(rff{Ej,)f{-Ej,). (28) 

The additional (21 + 1) comes from the degeneracy of the I 
states when the initial sum over the m-quantum number was 
taken. 

Since we are using 'trap units' the length and energy will 
be measured in harmonic oscillator length fl/,o = [h/mno)]^^^ 
and fiCL) respectively, and the temperature is in units of T = 
Hoj/kB- Using the LDA prescription we obtain the Fermi 
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FIG. 1. Spatial profiles of the (a,c) pairing field [A(r)] and 
the (b,d) the relative fluctuations in the pairing field [/(r) = 
-^((5"'"(r)(5(r))/A(r)]. In each of the plots three temperatures are con- 
sidered: T = O-OITf (solid blue curve), O-OSTp (dotted green curve) 
and O.lTp (dashed red curve) and the interaction strength is fixed at 
U = -5. For (a,b) = and for (c,d) P = 0.5. 

Energy Ep - (3Ny^^fia> and thus the Fermi temperature 
Tp - {3Ny^^fta)/kB- We can also calculate the Thomas-Fermi 
radius Rjp - {lANy^^Oho- The calculations were performed 
for = 10^ particles, at an interaction strength of U = -5 
(l/fepflj = -0.22) and a cut-ofF of E^. = ISOSw was chosen 
such that an increase in E^ does not change the converged re- 
sult. We have found that neglecting the Hartree terms does not 
qualitatively change the results presented here. Additionaly 
we only focus on the pairing field fluctuations, Eqs. (17) and 
(26), having found that the relative fluctuations in the densities 
Eqs. (18, 19, 27, 28) do not give any additional information 
concerning the validity of the meanfield solution. 

For a balanced system {P = {N^ - Ni)/N - 0) we regain the 
standard result for the pairing amplitude [Fig. 1(a)], i.e. away 
from the center of the trap A(r) decreases monotonically. The 
corresponding relative fluctuations in the pairing field, defined 
as/(r)= V<^V)^W)/A(r), are plotted in Fig. 1(b). This plot 
shows that the relative fluctuations in the paring field, increase 
with temperature and become significant only when A(r) — > 0. 
Since the Bogoliubov-de Gennes equation is valid for f{r) <*: 
1 we can ascertain that the meanfield approximation used to 
obtain Fig. 1(a) is valid, except when A(r) 0. 

For finite polaraization {P - 0.5) the FFLO state manifests 
itself as an oscillation in the meanfield pairing amplitude, for 
r/Rjp > 0.5 [Fig. 1(c)]. This state emerges as the temperature 
is decreased. As for the balanced case we can investigate the 
validity of the meanfield approach through the determination 
of the relative fluctuations in the pairing amplitude [Fig. 1(d)]. 
Importantly we find that in the region the FFLO state emerges 
the meanfield approximation breaks down, i.e. f{r) > 1, for 
T = 0.01 Tp and 0.05rF. However, we do note, as in the case 
of the balanced gas, that as the temperature is decreased the 
relative fluctuations in the pairing amplitude are reduced. As 
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FIG. 2. Spatial profiles of the (a) pairing field [A(r)] and the (b) 
relative fluctuations in the pairing field [/(r) = ^/{SHr)6(r))/A(r)]. 
In each of the plots three temperatures are considered: T = 0.0 ITp 
(solid blue curve), O.OSTf (dotted green curve) and O.lTp (dashed 
red curve) with coiTesponding polarizations of P = 0.895, 0.874 and 
0.841. The interaction strength is fixed at = -5 



such, in this case, we expect at extremely low temperatures 
(r < 0.01 Tp) that the meanfield predictions of FFLO states is 
valid. 

Ref. [22] considered the case where at a critical polariza- 
tion the FFLO states extend to the core of the superfluid for 
zero temperature. We have performed similar calculations at 
finite temperatures {T = 0.01 Tp, 0.05rp and 0.1 Tp) in the 
strongly interacting regime of U - -5 to examine a polarized 
core, [Fig. 2]. We found for polarizations of P = 0.895 ( 
0.874) that at T = 0.01 Tp (0.05rp) the FFLO state extends 
to the core of the superfluid [Fig. 2(a)]. Additionally, for 
r = O.lTp and P = 0.841 there is no FFLO state [red dotted 
curve in Fig. 2(a)], as calculated from the meanfield analysis. 
This is in stark contrast to the zero temperature calculation 
[22], where an FFLO state extending to the core is predicted. 
For T - O.lTp {P - 0.841) the relative fluctuations about the 
meanfield (/(r)) are of order one throughout the superfluid 
[Fig. 2(b)]. As such it is difficult to justify the validity of 
this meanfield solution. For the temperatures T = 0.05rp and 
0.01 Tp, where the meanfield FFLO state exists, the average of 
f{r) over the range r < O.lSRjp is greater than 1. As in the 
case of the edge FFLO states, shown in Fig. 1, the relative 
fluctuations about the meanfield decrease as the temperature 
is reduced. Specifically, for the parameters considered calcu- 
lations show that the meanfield analysis only becomes valid 
for T < 0.001 Tp, an experimentally challenging regime. 

Overall we have provided a general formalism for deter- 
mining the validity of meanfield Bogoliubov-de Gennes solu- 
tions based on the retrospective determination of the fluctua- 
tions about the meanfield, f{r). We have used this method to 
show specifically, the validity of the FFLO state in a spher- 
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ically symmetric ultracold gas. The analysis of the mean- 
field approximation at finite temperatures and zero polariza- 
tion shows the fluctuations to be small and confined to the 
boundary of the superfluid, i.e. as A(r) — > 0, indicating the 
meanfield approximation to be valid. However for a non-zero 
polarization we have seen /(r) at the boundary become more 
pronounced and the relative fluctuations become significant, 
casting doubt on the meanfield results. When we examine the 
FFLO states at the core of the superfluid, the relative fluctu- 
ations extend with the modulating pairing field and are sig- 
nificant for temperatures lower than T - 0.01 Tp. Only when 
we go several orders of magnitude below the current exper- 
imental temperatures can one be convinced of the meanfield 
results for finite polarization. In this study we have looked at 
a spherically symmetric trap, however one could extend this 
to asymmetric systems to consider the role of geometry and 
the validity of the meanfield approximation. 
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